module rad_funcs_mod

  use const_mod

  implicit none

  private

  public planck
  public integrated_planck

  real(r8), parameter :: c        = 2.997925d+10   ! Speed of light in vaccum (cm s-1)
  real(r8), parameter :: h_plk    = 6.62607015d-34 ! Planck's constant (J s)
  real(r8), parameter :: k_bol    = 1.380649d-23   ! Boltzmann's constant (J K-1)
  real(r8), parameter :: sig_sb   = 5.67032d-8     ! Stefan-Boltzmann constant (J m-2 s K-4)
  real(r8), parameter :: c2       = c**2
  real(r8), parameter :: hc       = h_plk * c
  real(r8), parameter :: hc2      = h_plk * c2
  real(r8), parameter :: hc_k     = hc / k_bol
  real(r8), parameter :: h_c2_2   = 2 * h_plk / c2
  real(r8), parameter :: k_h      = k_bol / h_plk

  ! Plank function parameters
  real(r8), parameter :: sigma_pi = sigma / pi
  real(r8), parameter :: conc     = 15.0d0 / pi**4
  real(r8), parameter :: vmax     = log(huge(1.0d0))
  real(r8), parameter :: eps      = radix(1.0d0) ** (1 - digits(1.0d0))

contains

  pure elemental real(r8) function planck(t, wn) result(res)

    ! J m-2 s-1 steradian-1 um-1

    real(r8), intent(in) :: t       ! Temperature (K)
    real(r8), intent(in) :: wn      ! Wavenumber (cm-1)

    res = 2 * hc2 * wn**5 / (exp(hc_k * wn / t) - 1.0)

  end function planck

  elemental real(r8) function integrated_planck(t, wn1, wn2) result(res)

    ! Integrating Planck's Function over given wavelength band to get the theoretical
    ! maximum amount of total (broadband) radiation that can be emitted by an object.
    !
    ! If over all possible wavelengths, we get the Stefan-Boltzmann law.

    !
    !  _ wf2                  _   _ 4  _ x2
    ! |                  2 h | k T |  |       x^3
    ! |      B(T) dwn =  --- | --- |  |     ------- dx
    !_|  wf1             c^2 |_ h _| _|  x1 e^x - 1
    !

    real(r8), intent(in) :: t       ! Temperature (K)
    real(r8), intent(in) :: wn1     ! Minimum wavenumber (cm-1)
    real(r8), intent(in) :: wn2     ! Maximum wavenumber (cm-1)

    real(r8), parameter :: coef = h_c2_2 * k_h**4
    real(r8) x1, x2
    integer n

    if (t < 1.0e-4) then
      res = 0
      return
    end if

    x1 = hc_k * wn1 / T
    x2 = hc_k * wn2 / T

  end function integrated_planck

  subroutine optical_depth(nl, ka, zi, tau)

    integer, intent(in) :: nl
    real(r8), intent(in) :: ka(nl)
    real(r8), intent(in) :: zi(nl+1)
    real(r8), intent(in) :: tau(nl+1)

  end subroutine optical_depth

end module rad_funcs_mod
